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Abstract. Monte Carlo simulation using the standard single-spin flip algorithm often fails to 
sample over the entire configuration space at low temperatures for frustrated spin systems. A 
typical example is a class of spin- ice type Ising models. In this case, the difficulty can be avoided 
by introducing a global-flip algorithm, the loop algorithm. Similar difficulty is encountered 
in 0(3) Heisenberg models in the presence of biquadratic interaction. The loop algorithm, 
however, is not straightforwardly applied to this case, since the system does not have a priori 
spin-anisotropy axis for constructing the loops. We propose an extension of the loop algorithm 
to the bilinear-biquadratic models. The efficiency is tested for three different ways to flip spins 
on a loop in Monte Carlo simulation. We show that the most efficient method depends on the 
strength of the biquadratic interaction. 



1. Introduction 

Monte Carlo (MC) simulation is a powerful tool for investigating thermodynamic properties of 
classical spin models pQ. The standard single-spin flip algorithm is widely used because it is 
simple and applicable to systems with any type of interactions. However, it is well known that it 
suffers from slow relaxation in many cases. For example, the single-spin-flip dynamics exhibits 
critical slowing down near a continuous transition temperature where the correlation length 
diverges. Such difficulty is avoided by using a nonlocal, global update, such as the Swendsen- 
Wang cluster algorithm for Ising systems [2] or its Wolff's extension to Heisenberg systems [3J- 
These cluster algorithms accelerate the relaxation by flipping many spins at once, i.e., updating 
the spin configuration drastically. 

Besides the critical slowing down, the single-spin flip algorithm often suffers from dynamical 
freezing at low temperature (T) when the ground state has macroscopic degeneracy under the 
influence of geometrical frustration. A typical example is an antiferromagnetic Ising model on 
a pyrochlore lattice. The pyrochlore lattice is a three-dimensional frustrated structure given 
by a corner-sharing network of tetrahedra, as shown in Fig. [Ha). When the antiferromagnetic 
exchange interaction is limited to nearest-neighbor sites, no long-range ordering occurs down 
to zero T and the ground state has macroscopic degeneracy 0]. The degenerate manifold is 
identified by a collection of local constraints enforcing two spins pointing up and two spins 
pointing down in every tetrahedron, as exemplified in Fig. HJa). This 'two-up two-down' 



constraint is called the ice rule because of an analogy to the constraint on positions of protons in 
hexagonal ice [5j[6]. Similar situation was recently discovered in the so-called spin-ice systems, 
in which the Ising-like spins point along the local (111) axes and interact with each other 
by ferromagnetic exchange interaction and dipole interaction [TJ [8]. Because the degenerate 
'ice-rule' configurations are separated by large energy barriers of the order of the dominant 
interaction scale J, the single-spin flip does not work at low T< J, The difficulty is avoided by 
introducing a global flip called the loop flip, in which one reverses all Ising spins on a specific 
closed loop passing through tetrahedra [HI [TUJ, [TTJ [T2J Q2] ; the loop is chosen so that the spins 
are up and down (inward and outward in the spin-ice problem) alternatively along the loop [see 
the hexagon in Fig. 0Ja) as an example]. This loop flip enables to transform an ice-rule state to 
another ice-rule state bypassing the energy barriers. 

The difficulty remains even when the Ising discreteness is relaxed and spins can fluctuate, as 
long as the ground-state manifold retains a multivalley structure with large energy barriers. Such 
situation is seen in variants of pyrochlore antiferromagnets, such as classical Heisenberg models 
in the presence of the single- ion easy-axis anisotropy |14|, [15] and the biquadratic interaction |16j . 
In contrast to the Ising case, however, it is nontrivial how to define the loop with alternating 
spins. Moreover, the loop flip procedure is not unique in the Heisenberg spin case because of 
the continuous degrees of freedom. Recently, the authors extended the loop algorithm to the 
Heisenberg spin systems with single-ion anisotropy by defining 'colors' of spins as black and 
white, which is a natural generalization of up and down in the antiferromagnetic Ising case, 
in terms of the projection of spins on the anisotropy axis [15]. We tested different ways to 
flip spins on a formed loop and showed that the efficiency strongly depends on the method. 
In the models with the biquadratic interaction, however, because of the 0(3) spin rotational 
symmetry, we have no explicit anisotropy axis to project spins on. Therefore, an extension of 
the loop algorithm to this class of models is not straightforward. 

In this paper, we develop an extension of the loop algorithm to classical antiferromagnetic 
Heisenberg spin systems with biquadratic interaction in which the spin-ice type manifold emerges 
at low T. We propose a way to define the projection axis for constructing loops and test three 
different ways to flip a formed loop. We apply the algorithm to a nearest-neighbor bilinear- 
biquadratic model, and compare the efficiency of the three methods. Interested readers are 
referred to Ref. |17] for further application of the present algorithm to a bond-disordered bilinear- 
biquadratic model. 

2. Extension of the loop algorithm to bilinear-biquadratic spin systems 

2.1. Model 

In this section, we extend the loop algorithm to classical antiferromagnetic Heisenberg models 
with biquadratic interactions. We start with a Hamiltonian of a simple form: 




where S, denotes a classical Heisenberg spin at site i on the pyrochlore lattice [Fig. Ufa)] (we 
take | Si | = 1) and b is the coupling constant of the biquadratic interaction. Here we consider 
b > which favors collinear spin configurations. We note that such 'ferro'-type biquadratic 
interaction originates in the spin-lattice coupling as well as quantum and thermal fluctuations. 
We consider the antiferromagnetic exchange interaction J > 0, and take the energy unit as 
J = 1. The sum runs over the nearest-neighbor bonds of the pyrochlore lattice. The following 
algorithm is applicable to more general spin-ice type models on other frustrated lattices with 
farther-neighbor or bond-dependent interactions, and site-dependent anisotropy. 




Figure 1. (a) The pyrochlore lattice composed of a three-dimensional network of corner-sharing 
tetrahedra. A 16-site cubic unit cell is shown. Spins are denoted by solid arrows. The broken 
arrow Q shows a common axis selected by b in the nematic phase. Black (filled) circles represent 
spins in the direction of Q, while white (open) circles represent spins in the opposite direction. 
The spin configuration is an example of the spin-ice type states in which the 'two-up two-down' 
local constraint is satisfied in every tetrahedron. The hexagon with a bold dashed line denotes 
one of the shortest loops on which a flip of all spins (black and white) transforms the ice-rule 
state to another ice-rule state, (b) Three different ways to flip black and white: (1) flip xyz, (2) 
flip parallel, and (3) rotate. See the text for details. 



When 6 = 0, the model given by Eq. ([T]) exhibits no long-range ordering down to T = 0, and 
the ground state is given by a collection of local constraints that enforces the sum of Si to be 
zero in every tetrahedron. Consequently, the ground-state manifold has continuous macroscopic 
degeneracy [18^ [19l I20| . For b > 0, the present model exhibits a weak first-order transition 
at T c ~ b to a nematic state in which spins spontaneously select a common axis Q without 
selecting their directions on it |16j . Hence, the ground state for b > is identified by a collection 
of spin-ice type local constraints: in every tetrahedron, two out of four spins are aligned parallel 
to each other and the other two are antiparallel to them — 'two-up two-down' configuration [see 
Fig.QJa)]. The ground-state manifold develops a multivalley structure whose minima correspond 
to different spin-ice-type configurations separated by large energy barriers of the order of b and 
J. 

2.2. Algorithm 

In the previous paper, the authors developed an extension of the loop algorithm to be applicable 
to a family of classical Heisenberg antiferromagnets with single-ion anisotropy, in which the 
single-spin flip algorithm suffers from slow relaxation due to the formation of the spin-ice type 
manifold [15] . In the extended loop algorithm, (i) we first project all the spins onto the anisotropy 
axis to assign black and white colors, (ii) next, construct a loop of alternating black and white, 
and (iii) flip all the spins on the loop. For the present bilinear-biquadratic model, a similar 
difficulty from slow relaxation is anticipated because the low-T state develops the spin-ice type 
degeneracy as mentioned above. However, it is not straightforward to apply the extended loop 
algorithm since the present model retains 0(3) spin rotational symmetry and does not have any 
explicit anisotropy axis to project the spins on. It is necessary to deduce the common axis Q 
selected by b for each MC sample. 

Here, we propose the following procedure to define the projection axis. We first pick up a 
set of iVx tetrahedra {7~ m } (m = iVr) randomly from the whole system. Starting from an 



initial guess «o [we take «o = (0,0, 1)], the normalized projection axis a is obtained iteratively 
by 



a n+1 oc sign(5j • a n )Si. (2) 

ie{T m } 

Here the sum is taken over all spins belonging to the selected tetrahedra {T m }, and n (= 
0, 1, • • • , n max — 1) is the index of the iteration. For larger iVr and ri max , the resultant a = a nmax 
gives a better approximation of Q. In practice, we take iVr = 24 and n max = 6 in the following 
calculations for the system size with 7V S = 16 x 8 3 spins. We confirm that the loop flip is 
efficiently performed for these conditions, as demonstrated below. Once a is defined in this 
way, the step (ii) and (iii) are done in the same way as in Ref. |15j . It should be noted that, to 
ensure the detailed balance, loops must be constructed avoiding the tetrahedra included in {T m } 
as well as defect tetrahedra in which the ice rule is violated: Otherwise, the loop flip becomes 
irreversible because the flip changes a. 

In the extended loop algorithm, the loop flip procedure to reverse all colors on a loop is 
not unique because spins can thermally fluctuate [15]. To choose an efficient method, careful 
consideration on the energy change is necessary. In Ref. [15], the authors tested two different 
ways: (1) flip xyz and (2) flip parallel. In flip xyz, all three Cartesian components of Si are 
reversed as Si — > —Si, while in flip parallel, only parallel components S^ are reversed as 

Si — > Si — 2(Si ■ a)6i [see Fig. Ifb)]. The previous study revealed that, in the case of the 
single-ion anisotropy, only the acceptance rate for flip parallel can become one (rejection free) 
in the limit of T — > 0; the acceptance rate for flip xyz converges to a finite value less than unity 
as T — > because of the effect of thermal fluctuations on the transverse component of spins |15j . 
In addition to these two updates, in this paper, we introduce another way to reverse colors, i.e., 
(3) rotate. In rotate, one translates every spin to the neighboring site on the loop simultaneously 
in the direction in which the loop was formed [see Fig. QJb)]. In the next section, we try the 
three different methods and demonstrate that the most efficient method depends on the model 
parameter b. 

3. Benchmark 

In this section, we apply the extended loop algorithm to the model given by Eq. ([T]). We 
demonstrate the efficiency of loop flips in MC simulations and compare the efficiency of the 
three methods. In the following, we show the MC results for the systems size with N B = 16 x 8 3 
spins under periodic boundary conditions. To retain the ergodicity, we use the loop flip together 
with the single-spin flip. One MC step consists of single-spin flips, followed by loop flips with 
either flip xyz, flip parallel, or rotate. In the single-spin flips, we randomly choose a new spin 
state on the unit sphere for each spin following a procedure proposed by Marsaglia [21] , 

At low T compared to b and J, spin configurations are enforced to satisfy the 'two-up two- 
down' ice rule, and the acceptance rate of the single-spin flip, P s i ng i e , is suppressed. This is 
demonstrated in Fig. [2l On the other hand, the probability that a closed loop is successfully 
formed, -Pi 00 p> steeply increases below the nematic transition temperature T c ~ b, indicating 
that almost all tetrahedra start to follow the ice rule below T c . [For b > J, the ice rule is weakly 
violated in the range of T > J even below T c , but gradually satisfied for T < J, as exemplified 
in Fig. [2]^d).] The acceptance rate of flips of a formed loop also increases below T c and remains 
finite as T — > 0; here, P xyz , -Pparaiieh and -Protate are the rate for flip xyz, flip parallel, and rotate, 
respectively. The total acceptance rate of the loop flip is given by the product as Pi 00 p x Pxyz, 
-Pioop x ^parallel) and P r otate x -^parallel for each method. Hence the acceptance rate of the loop 
flip sharply increases at T < T c , compensating the decrease of P S mgic- These T dependences are 
qualitatively the same for the wide range of b, as shown in Fig. [2l 




Figure 2. Temperature dependences of the acceptance rates of the single-spin flip (-P S mgie); the 
probability of formation of closed loops (P\ 00 p), the acceptance rates of flip of a formed loop 
by flip xyz (P xyz ), flip parallel (-P pa raiiei)> and rotate (-Protato)- The data are calculated at (a) 
b = 0.05, (b) 0.2, (c) 0.6, and (d) 1.5. The vertical broken lines denote the nematic transition 
temperature T c estimated by the peak position in the specific heat (not shown). 




As demonstrated in the previous study for the models with single- ion anisotropy [TB], the 
efficiency of the loop flip at low T depends on the method. Furthermore, in the present case 
with the biquadratic interaction b, the efficiency at low T strongly depends on b. The result is 
presented in Fig. [3j For small b < 0.1, flip parallel is most efficient, while it is taken over by 
rotate in the intermediate range of 0.1 < b < 0.5, and finally by flip xyz for b > 0.5. 

The difference of the efficiency is understood by the following consideration on the energy 
change by the flips. Considering a given state at a finite T well below T c , its energy measured 
from the ground-state energy is given by E = Ej + Ef,, where Ej and E^ are the energies 
corresponding to the first and second terms in Eq. (TjQ), respectively. Both Ej and E^ are of 
the order of T at low T. The three loop flips change the two contributions in different ways. 
The flip xyz changes Ej by a certain fraction AEj but conserves Ef,. Meanwhile, flip parallel 
and rotate change both of Ej and Ef, by AEj and AE},, respectively. First we consider the 
large b limit, where E ~ Ei, oc T and Ej/Eb oc b . For flip xyz in which AE^ = 0, we obtain 
AE/T = AEj/T oc 6" 1 — > as b — > oo. Since the acceptance rate is given by exp(— AE/T), 
this consideration gives lim6^. +00 P xyz — > 1, that is, flip xyz becomes rejection free as b — >■ +oo. 



This is consistent with the behavior in Fig. [3l On the contrary, flip parallel and rotate cannot 
become rejection free for b — > oo because they change Ef, by 0(T); P pa raiiei and -Protate converge 
to finite values less than unity at b — >■ oo, respectively. Note that rotate does not change a half 
of the nearest-neighbor bond energies. This may account for why P ro t a te > -fparaiiei at b — > +00. 

For smaller b, spins deviate from the common axis Q with larger angle. Considering that 
Q is set by Ei ~ T and Et, is proportional to bQ 2 (6 is a typical deviation angle), we obtain 
9 = Oi^fTjb) for b < J [22]. Since flip parallel changes Ej and £ 6 by 0(Jfl 4 ) and 0{b9 2 ), 
respectively [TS], we obtain AE/T = AE^/T = O(l) at low T <^ T c for parallel. This 
indicates that -P pa raiiei does not vanish even at b —> in the nematic phase. Meanwhile, the 
other two methods change the energy as AE/T = AEj/T = 0(J6 2 /T) = 0(J/b) at low T and 
b <C J [IS]. This suggests that their acceptance rates vanish as b — > in contrast to flip parallel. 
Therefore, flip parallel becomes most efficient at b — > 0. These considerations are consistent 
with the numerical results shown in Fig. 3. In the intermediate regime, i.e, 0.1 < b < 0.5, rotate 
is superior to the other two, presumably because of a remnant of the advantage of rotate over 
flip parallel at b — > 00. 

4. Summary 

In this paper, we have extended the loop algorithm to the classical antiferromagnetic Heisenberg 
spin models with biquadratic interaction which have spin-ice type ground-state degeneracy. The 
efficiency of the extended loop algorithm has been demonstrated in Monte Carlo simulations. We 
have examined three different ways of loop flips, flip xyz, flip parallel, and rotate, and compared 
their efficiency. We have shown that the most efficient method depends on the strength of 
the biquadratic interaction b. This b dependence has been explained by considering effects of 
thermal fluctuations on the energy changes by the loop flips. 

This work was supported by Grant-in- Aids (No. 19052008), Global COE Program "the 
Physical Sciences Frontier", and the Next Generation Super Computing Project, Nanoscience 
Program, MEXT, Japan. 
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